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Abstract 


Nuclear energy is emerging as a major alternative to the rapidly depleting conventional sources of 
energy. Safety is one of the major concerns of the nuclear power industry. The pressure vessel 
containing the nuclear reactor assembly has a number of inlet and outlet no22les welded to it. It is 
important to demonstrate the structural integrity of the containment vessel to ensure a safe 
operation. One of the most important aspects in this regard is the consideration of the surface 
cracks in the crotch region of the nozzle-vessel intersection, called nozzle-comer cracks. 

The crack characterizing parameter, stress intensity factor, is used to estimate the severity 
of a crack, within the framework of Linear Elastic Fracture Mechanics. Due to the complexity of 
geometry, the stress intensity factor for nozzle-comer cracks can only be evaluated using either 
approximate methods, or more accurately, by using numerical techniques such as the Finite 
Element method. 

In the present work, the Finite Element method has been used to evaluate the stress 
intensity factor variation along the profile of a nozzle-comer crack, modeled as a quarter ellipse, in 
a given nozzle-vessel intersection. Two types of load conditions have been considered: A simple 
internal pressure load, representing the normal operational condition; and a pressurized thermal 
shock, representing a faulted or accidental condition. 

Firstly, for the normal operational condition, a three-dimensional static stress analysis of 
the nozzle-vessel intersection without crack has been performed to demonstrate the validity of the 
partial modeling of the structure and the boundary conditions used. The comer crack is then 
included in the model and the stress analysis has been repeated, with singular elements near the 
crack front. Stress intensity factor has been evaluated from the displacement field near the crack 
profile, as obtained from the stress analysis. The results have been validated by solving two test 
problems and also by comparison with an empirical formula. Parametric study has been done for 
different crack and geometrical parameters. 

Secondly, for the faulted condition, a given crack configuration, with crack depth equal to 
half the vessel thickness, has been studied. A transient thermal stress analysis has been performed 
and the stress intensity factors have been evaluated at different points in the time domain, and are 
compared with the critical stress intensity factor for crack initiation. The particular half-thickness 
crack considered here has been estimated to be critical. 
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Further, as an attempt towards the estimation of stable growth of nozzle-comer cracks, a 
three-dimensional formulation of the cohesive-zone method for crack growth analysis has been 
presented. The method has been applied for the crack growth analysis of a Double Cantilever 
Beam Specimen and the results have been found to match qualitatively well with the available 
literature. The method is then applied for the estimation of growth of a given nozzle-comer crack 
under normal operational load and the pattern of crack growth has been found to be consistent with 
the SIF results obtained. 
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Chapter 1 


INTRODUCTION 


1.1 General Introduction 

Pressure vessels are leak-proof containers supposed to sustain high temperature and pressure, 
usually under adverse environmental conditions, which typically include corrosion, neutron 
bombardment, hydrogen embrittlement, etc. The ever-increasing use of vessels for storage, 
industrial processing, power generation, chemical industries, space and nuclear applications has 
given special emphasis to analytical and experimental methods for determining their operating 
stresses. High pressures, extremes of temperature, severity of functional performance 
requirements, eind above all, safety requirements, pose exacting design problems. One of the key 
aspects of design is the identification of the most likely modes of damage or failure. 

Especially, when it comes to applications like nuclear power generation, the requirements 
of safety increase manifolds. With the rapid depletion of the conventional sources of energy, 
nuclear energy is emerging as a major alternative. Owing to the abundant thorium reserves, India 
has a huge potential for the nuclear power industry. This makes the integrity of containment 
systems of nuclear reactors, an important field of study. One of the most significant modes of 
failure of nuclear reactor pressure vessels is fracture. Fracture could be initiated at existing defects 
or at cracks grown to a critical size by fatigue. 

Pressure vessels usually have a number of inlet and outlet nozzles attached. Under the 
action of internal pressure, there exists a high stress concentration at the inner comer of the nozzle- 
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vessel intersection. Besides, due to the difficulty of welding, it is possible to produce initial flaws 
at the intersection. This makes the nozzle crotch comer region susceptible to damage through 
fatigue and subsequent failure through fracture. Hence, the fatigue and fi:acture of nozzle-comer 
cracks is one of the most important considerations in the assessment of integrity of pressure 
vessels. 

There are two types of fracture that can be observed in a single crystal; brittle fracture and 
ductile fracture. A brittle fracture takes place without substantial plastic flow and occurs when the 
normal stress on one of the principal planes of the crystal reaches the critical value. When large 
plastic deformation occurs, consisting of sliding along crystallographic planes, prior to fracture, the 
failure is called a ductile fracture. These two types of fracture are also characteristic of 
polycrystalline materials. The brittle type exhibits little deformation and has a flat fracture surface 
essentially perpendicular to the direction of maximum stress; whereas, the ductile type exhibits 
gross deformation and the fracture surface has the familiar cup-and-cone form with shear-lips. 

The thicker a component, the greater is the tendency for brittle, rather than ductile fracture 
to occur. Hence, linear elastic fracture mechanics (LEFM) theory can be used to assess the severity 
of cracks in thick pressure vessels, such as those employed in light water reactors (LWR). To 
ensure the structural integrity of a pressure vessel, it must be shown that sufficient safety margins 
exit between crack sizes which can be safely detected and those crack sizes which can lead to 
failure, under operational, as well as faulted conditions. 

The crack characterizing parameter, stress intensity factor (SIF), proposed by Irwin in 
1957, is used for this purpose within the framework of LEFM. The stress intensity factor is an 
important quantity in both fatigue and brittle fracture situations. In the former, the SIF range can be 
related to fatigue crack growth and in the latter, SIF can be compared with the material fracture 
toughness (critical stress intensity factor) to asses the likelihood of brittle fracture. Although the 
evaluation of SIF is fairly simple for cracks in the cylindrical part of the vessel, the situation is 
much more complicated for nozzle-comer cracks. Because of the complexity of the geometry, it is 
not possible to derive analytical solution. It can be estimated only by an approximate method. The 
approximate methods that have been used so far to estimate the value of SIF can be roughly 
categorized in to the following three kinds: 
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1) Approximate analysis: Among these, Kobayashi’s approximate analysis, Besuner’s 
influence function method and Mohamed’s approximate method are the most well 
known. However for all these methods, the stress distribution at the uncracked geometry 
must be known beforehand. 

2) Experimental methods: These involve mainly the fatigue crack growth method, photo- 
elasticity method etc. The advantage of these methods is that their results are reliable as 
they involve the effect of all the factors on the SIF. But tiie shortcomings are the effort 
and cost in model making and testing. 

3) Finite element method (FEM): The main advantage of this method is its general purpose 
nature and that it can be used to study the cracks of any shape under any kind of loading. 
Its results are relatively reliable. However it is costly in computation. 

Since, none of these methods can fetch the value of SIF instantly, it is desirable to apply one or all 
of these, and compile SIF data for a number of no 2 zle-vessel geometries in a given range, and a 
number of crack shapes and sizes. 

Besides the consideration of the normal operational load condition of the Reactor Pressure 
Vessel (constant internal pressure load), it is also important to consider the accidental conditions. 
Under normal circumstances, the temperature gradients in the nozzle-vessel intersection are such 
that they produce compression on the inner surface, and hence the thermal stresses do not cause 
any adverse effect on a nozzle-comer crack. Therefore, in the normal operational condition, it is 
sufficient to consider only the internal pressure load to study the behaviour of nozzle-comer 
cracks. 

But in a special situation, commonly referred to as the “thermal shock” problem, which 
arises during a typical accidental condition due to a sudden loss of the Reactor coolant, the thermal 
stresses can prove fatal for a nozzle-comer crack. Thermal shock involves the possibility of 
fracture of a nuclear reactor pressure vessel during a loss of coolant accident (LOCA) under 
circumstances that are most likely to occur in a pressurized water reactor (PWR) plant. During a 
LOCA, cold water is splashed on the inner surface of the vessel by the emergency core cooling 
system (ECCS). This produces thermal stresses which can cause a nozzle-comer crack, which is 
otherwise safe, to become critical. For this reason, the thermal shock problem is of profound 
concern to the nuclear power industry. 
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1.2 Literature Review 


Three dimensional finite element methods have primarily been employed to determine stress 
intensity factors for a crack at the crotch of a nozzle in a pressure vessel. The first SIF solution for 
such cracks was published by Rashid and Gilman [1], who based their finite element solution on 
strain energy release rate. 

Shah and Kobayashi [2] gave the solution of an embedded elliptical crack exposed to 
arbitrary normal pressure. 

Hellen and Dowling [3] gave the SIF solutions around the profiles of four crack 
configurations at the inner crotch of a typical LWR no 22 ;le-vessel intersection. Three dimensional 
finite element techniques were used, incorporating quadratic isoparametric elements with special 
singularity elements around the crack profile. Direct substitution in the displacement field solution 
near the crack tip, and virtual crack extension techniques were used for the calculation of SIF. 

Schmitt [4] used a numerical technique based on the Westergaard stress function to 
calculate the value of SIF at mid points of straight-front cracks of the nozzle of a vessel. 

Schmitt [5] analysed a nozzle-comer crack by FEM using crack tip singularity elements 
and compared the SIF values along the crack fi-ont with the crack initiation fi:acture toughness 
{K,^) and the crack arrest fracture toughness ( ) to demonstrate the structural integrity. 

Broekhoven [6] used so-called semi-analytical procedure to find SIF for the cracks at the 
comer of a nozzle-plate intersection. The technique is a modification of the solution of an 
embedded elliptical crack by Shah and Kobayashi [2]. 

Bessuner et al. [7] employed the influence function (IF) method to calculate average SIF 
for cracks of various sizes and locations at the inside of a nozzle attached to a vessel. The method 
involves the use of an influence function which is independent of the loading but depends only on 
the crack shape, nozzle geometry and boundary conditions. 

Kobayashi et al. [8] has obtained SIF values for nozzle-comer cracks by employing the 
results of a quarter-elliptical crack at a quarter space exposed to different loading conditions. These 
results were obtained using the so-called alternating technique. 

Mohamed and Schroeder [9] used the existing solutions for comer cracks at the 
intersections of nozzles with plates or vessels to develop an empirical relationship for SIF of such 
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cracks. They also proposed a general method to predict SIF values of comer cracks using stress 
concentration for pressure loading only. 

W. Schmitt et al. [10] gave the SIF solutions for two different nozzle geometries and 
different crack sizes for pressure and thermal loading using three dimensional elastic finite element 
models. They also proposed a procedure to estimate the maximum SIF for arbitrary crack size and 
loading conditions. 

Guozhong and Qichao [11] simplified the iso-stress lines at a nozzle-comer into straight 
lines at an angle of 45 “ with the walls of the cylinder and nozzle, and based on this, derived an SIF 
solution for an arbitrary point on the crack front of a nozzle-comer crack by means of an 
approximate analysis. Following from this, they proposed approximate expressions for maximum 
and average SIF. 

Baisong et al. [12] proposed an improved 3-D collapsed isoparametric singular element 
with quarter-points for crack analysis using FEM. They also analysed a comer crack in the interior 
wall of a nozzle of a Reactor pressure vessel using the proposed element, and calculated the SIF 
using the William’s formula for displacement field near a stationary crack, as well as the least- 
squares linear fitting method. 

The only known experimental results for nozzle-comer cracks are those of Derby [13], 
Broekhoven [14], Smith et al. [15] and Ruiz [16]. Derby [13] conducted burst tests on two groups 
of epoxy models of intersections of nozzles with vessels. The average SIF at burst was taken to be 
equal to the fracture toughness of the material. Smaller sizes of the models were also investigated 
experimentally by Smith et al, [15] using the frozen-stress photoelastic technique, which permits 
the evaluation of variations in SIF values along the crack fronts. 

Broekhoven [14] correlated the propagation rates obtained for nozzle-comer cracks with 
those of standard specimens of the same material. He obtained minor variation of SIF values on the 
crack front. 

Ruiz [16] also used the frozen-stress photoelastic technique to determine the average SIF 
for straight-front cracks at the comer of several models of tees (intersection of two cylinders). 

Smith et al. [17] used a method which is a combination of frozen stress photoelasticity and 
a computerized least-squares data analysis for measuring SIF distributions at nozzle-comer cracks. 
The authors have emphasized on the importance of using actual flaw shapes in analysis. 
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1 .3 Overview of the Present Work 


The problem considered in the present work is given for study by the Pressure Vessels Section of 
the Light Water Reactor Division (LWRD) of Bhabha Atomic Research Centre (B.A.R.C.), 
Mumbai. The problem is concerned with the assessment of structural integrity of a typical, thick, 
high pressure vessel having a nozzle-comer crack. Although, design codes, such as the ASME 
Boiler and Pressure Vessel code [18], describe methods to roughly estimate SIF for cracks in 
pressure vessels, they are often over conservative and do not reveal the true picture; and 
experiments, as stated previously, are costly and more time consuming. In lieu of the above, 
numerical methods like FEM offer a viable and reliable solution. 

The present work can be subdivided in to three main parts, as detailed below, along with a 
brief description of each part; 

I. Calculation of Stress Intensity Factors for nozzle-comer cracks under normal operational 
load condition: 

1) Stress analysis of the uncracked nozzle-vessel intersection and its validation. 

2) Stress analysis of the intersection with comer crack. 

3) Evaluation of SIF along the crack front using the results of stress analysis, and its 
validation. 

4) Parametric study of SIF of nozzle-comer cracks for different geometrical and crack 
parameters. 

II. Study of the behaviour of a given nozzle-comer crack during a LOCA: 

1) Transient thermal-stress analysis of nozzle-vessel intersection with comer crack. 

2) Evaluation of SIF along the crack front at different time points. 

3) Plotting of SIF at different points of the crack front as functions of time and their 
comparison with material fracture toughness for crack initiation. 
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III. Prediction of the amount of stable crack growth using cohesive-zone method; 

1) Incorporation of the special cohesive elements in to the Finite Element (FE) model 
and inclusion of the corresponding FE equations in to the general purpose FE code 
ABAQUS, through a user subroutine written in FORTRAN. 


1 .4 Organization of the Thesis 


The thesis consists of seven chapters in all. The organization of the thesis is as follows: 

1 ) Chapter 2 gives the FE formulations of 3D stress and heat transfer analyses. 

2) Chapter 3 presents the results and validation of stress analysis of nozzle-vessel intersection 
without crack. 

3) Chapter 4 deals with the effect of the normal operational load on nozzle-comer cracks. The 
Chapter presents the results of stress analysis of the nozzle-vessel intersection with a 
comer crack. The SIF is then evaluated for the comer cracks using the method described in 
the Chapter and the results are verified. This is followed by the parametric study of SIF. 

4) Chapter 5 is dedicated to the study of effect of a LOCA on a nozzle-comer crack. 

5) Chapter 6 gives the formulation of a 3D cohesive zone element and presents the results of 
crack growth analysis. 

6) Chapter 7 concludes the work and suggests the scope for future work. 
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Chapter 2 


FE FORMULATIONS FOR STRESS AND HEAT 
TRANSFER ANALYSES 


2.1 Introduction 

The study of a physical phenomenon usually involves two major tasks: Mathematical modeling of 
the physical process resulting in algebraic, differential or integral equations and solution of the 
resulting equation. 

While the derivation of the governing equations for most problems is not unduly difficult, 
their solution by exact methods of analysis is a difficult task. In such cases, approximate methods 
of analysis provide alternative means of finding solutions. 

The Finite Element (FE) method is one of the most popular methods of numerical analysis, 
in which, a given domain is represented as a collection of simple sub-domains, called finite 
elements, so that, it is possible to systematically construct the approximation functions over each 
element, needed in a variational or weighted-residual approximation of the solution of a problem. 

In this chapter, the FE formulations for a 3D static stress analysis, and a 3D transient heat 
transfer analysis have been presented. 


2.2 FE Formulation for 3D Static Stress Analysis 

Consider an arbitrary three-dimensional solid body Q with volume V and boundary F as shown 
in Figure 2.1. 
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Figure 2.1 Domain and Boundary Conditions for Stress Analysis. 



Figure 2.2 Domain and Boundary Conditions for Heat Transfer Analysis. 
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With reference to the Figure, represents the boundary on which displacement is 
specified. Traction t is specified on the boundary F, and represents the concentrated forces 
acting on the body. 

The principle of virtual work states that for a virtual perturbation in the displacement field 
at equilibrium, the virtual change in the internal strain energy {8U) must be balanced by an 
identical change in the external work {SW ) due to the applied loads, that is: 

SU^W, ( 2 . 1 ) 


or. 




( 2 . 2 ) 


where, SU^ and are the contributions of the finite element e to the virtual change in strain 

energy and external work respectively. 

Now, 


5U,= \{5sY{a]dV^, (2.3) 

where, represents the volume of the element e , 

{ff} = Strain vector = [f, s, , 

{cr}= Stress vector = [o-, cr^ cr, cr^ <x^ <jJ[. 

Assuming an isotropic linear elastic material behaviour. Equation 2.3 may be expanded as: 

(2.4) 

I'e 

where, [f>] is the material property matrix given by: 
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E is the Young’s modulus of the material and v is the Poisson’s ratio; 


and { 5 ,^}= Thermal strain vector = a a 0 0 Oj , (2.5) 


where, AT = T-T^, 

a = Coefficient of thermal expansion of the material, 

where, T = Current temperature at the point in question, 

= Reference (strain-free) temperature. 


The strains are related to the displacements as: 


du 

dv 

dw 

TSu dv^ 

dx 


dz 



dv dw^ 
dz dy^ 


du ^ frw'’ 
dz dx j 


( 2 . 6 ) 


where, « , v and w represent the displacements of any arbitrary point within the element, in jc , y 
and z directions respectively, and are interpolated from the nodal displacements through the 
element shape functions as: 




[W 


(2.7) 
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where, [A^] is the element shape function matrix and is the element nodal displacement 
vector. 

From Equations 2.6 and 2.7, 


{s}=[b]{u,}, ( 2 . 8 ) 

where, [5] is the strain-displacement matrix based on the derivatives of the element shape 
functions. 

Combining Equations 2.4 and 2.8, and noting that {u^] does not vary over the volume of the 
element: 




(2.9) 


Now, 


■w. ={«»,)" pYk')>r; 


r,‘ 


( 2 . 10 ) 


where, F" is the area of the element that falls on F^, 

= Traction vector on F‘ = [/' r * tl , 

(p/]= Concentrated nodal force vector. 

Finally, substituting Equations 2.9 and 2.10 into Equation 2.2 and writing the Equation on element 
basis to give: 

£( {((»,}'■ ^Bf[D][B\iV,{«,}-lSuJ JisrWkK )-E( \Su,Y J[jvr + ) 

« K. K, e rf 
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Now, since is arbitrary, Equation 2. 1 1 reduces to: 


SKIKl'Ik'Kik'KEk*). (2-12) 

e e e e 


where, \k^ ] = J[5]^ = Element stiffness matrix, 

I'e 

{^e }“ Jt^r - Element traction load vector, 

r,‘ 

{?'/*}= = Element thermal load vector. 


The element matrices and load vectors of all the elements in the domain are assembled together to 
get the global stiffhess matrix [iiT] and the global load vector {f} such that: 

Mk)=W. 

where, } is the global nodal displacement vector. 


2.3 FE Formulation for 3D Transient Heat Transfer Anaiysis 

The governing differential equation for transient 3D heat transfer in a solid is: 

pC^ + {LY{q]^q, (2.13) 

at 

where, p = Density of the material, 

C = Specific heat of the material. 
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T = r(x,>’,z,/)=Temperature, 


t =Time. 


{!} = 


dx 

dy 


dz) 


= Vector operator, 


= Heat flux vector, 

q = Heat generation rate per unit volume. 


Fourier’s heat conduction law yields: 


where. 


{q} = -[D]{L}T, 


( 2 . 14 ) 


[Z)] = 


K 0 0 
0 0 
0 0 k. 


= Conductivity matrix. 


where, k^,ky and k^ are thermal conductivities in x, y and z directions respectively. 
Combining Equations 2.13 and 2.14, 


pC^j^-{LY([D]{L}TM. 


( 2 . 15 ) 


Three types of boundary conditions are considered (with reference to Figure 2.2); 

1 ) Specified temperatures over surface S , : 

T = T\ (2.16) 

where, T‘ is the specified temperature. 
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2) Specified heat flux over surface : 



(2.17) 

where, {«} is the unit outward vector normal to the surface and q’ 

is the specified heat 

flux. 


3) Specified convection over surface : 


{«}’■{»)= 

(2.18) 

where, hj- is the convective film coefficient, Tg is the bulk temperature of the adjacent 

medium and is the surface temperature. 


Combining Equations 2.14 with Equations 2.17 and 2.18, 


{nY[D]{L}T = q\ on S,, 

(2.19) 

{nY[D]{L}T^hjiTg-T), on S,. 

(2.20) 


Pre-multiplying Equation 2. 1 5 by a virtual change in temperature, integrating over the volume ( ) 
of the element e, and combining with Equations 2.19 and 2.20, with some manipulation yields the 
variational statement: 

+ (5r)([D]{i}r)lt/F, = + jsrh^{T, -T)dS, + j^qdV^ , (2.21) 

) Si S3 Ke 

where, dT = Sr{x, y, z, t) = Virtual change in temperature. 

Now, the dependencies of T on space and time are separated as: 

r = {vy'{7;}, (2.22) 
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where. 


{n} = {//(jc,;;, z)} = Element shape function vector, 
{T^ } = (r )} = Element nodal temperature vector. 


Thus, 




(2.23) 


sr^{^J[N}, 

(2.24) 


{i)r=[s]{2;}, 

(2.25) 

where. 




Now, Equation 2.21 can be combined with Equations 2.22 - 2.25 to yield: 

‘ ye ‘ S2 

(2.26) 

+ \\^,Y{N]h,{f,-{NY{T,'^dS,*l\^,Y[N}qdV, ). 

S3 Ve 


p is assumed to be constant over the volume of the element, while C and q may vary over the 
element. Now, since {<5r^} is arbitrary, Equation 2.26 reduces to: 

K [C.]{7;}+(kl+kl)fc} )=2( fc^l+tel+b) ). (2.27) 

e e 


where, [C'J = /? - Element specific heat (thermal damping) matrix, 

ye 
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[^c]~ jl^] Element diffusion conductivity matrix, 

‘i- 


[^‘ ] = = Element convection surface conductivity matrix, 


= Element heat flux vector, 


{Se } == l^B^f {^]^3 - Element convection surface heat flow vector, 

i'i 


{2/}= = Element heat generation load vector. 

Ve 


The element matrices and load vectors are assembled together for all the elements in the domain to 
get the global thermal damping matrix ([c]), the global conductivity matrix ([^]) and the global 
load vector ( {q}) such that; 


[c]{rJ+Mfc)={e}. 


where, } is the global nodal temperature vector. 


2.4 Closure 

The Finite Element formulations for three-dimensional static stress analysis and transient heat 
transfer analysis have been presented. In the present work, static stress analysis of the nozzle- 
vessel intersection has been performed for the normal operational load condition; and a sequential 
transient thermal-stress analysis has been performed for the faulted condition (LOCA). 
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Chapter 3 


STRESS ANALYSIS OF UNCRACKED NOZZLE-VESSEL 
INTERSECTION 


3.1 introduction 

The aim of this chapter is to identify the partial section of the structure in the vicinity of the 
nozzle-vessel junction, which is suitable and sufficient for the purpose of stress analysis, and also 
to validate the overall procedure adopted for the stress analysis. The normal operational load 
condition (simple internal pressure load) has been chosen for the purpose. Stress analysis of a 
partial section of the structure (without crack) has been performed and the methodology has been 
validated by comparison with the work of Moini and Mitchell [19]. 


3.2 Geometry 

The basic geometry of the nozzle-vessel intersection considered in the present work, as given by 
the Pressure Vessels Section of LWRD, B.A.R.C., is described below: 

Figure 3.1 shows the nozzle-vessel intersection, sectioned by the plane which contains the 
axes of both, the vessel and the nozzle. The main dimensions characterizing the intersection are 
also shown in the figure. 
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Figure 3.1 Geometry of the nozzle-vessel intersection. 


These dimensions are given as: 

1 ) Inner radius of the vessel, R, = 700 mm. 

2) Outer radius of the vessel, R 2 = 1000 mm. 

3) Inner radius of the nozzle, r,= \07 mm. 

4) Outer radius of the nozzle, = 165 mm. 

5) Radius of curvature of the inner fillet, r, = 50 mm. 

6) Radius of curvature of the outer fillet, r^= 100 mm. 


3.3 Material 

Nuclear grade low alloy ferritic steels are commonly used for the fabrication of nuclear reactor 
pressure vessels. The material considered in the present work is the ASME SA 508-Class 3 
pressure vessel steel. The chemical composition of the steel [20] is given in Table 3.1. 
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Element 

Wt. % 

Element 

Wt. % 

Carbon 

0.18 

Silicon 

0.03 

Manganese 

1.48 


6.004 

ESi&luOHli 

0.003 

Nickel 

0.91 

Chromium 

0.21 


0.54 

Aluminium 

0.003 


0.045 

Vanadium 

0.003 

Cobalt 

6.004 


Table 3.1 Chemical composition of SA 508-Class 3 steeL 


The relevant mechanical and thermal properties of the material have been taken from Sattari-Far 
[21], where, the material properties for SA 533 B steel is given, which is similar to SA 508-Class 3 
steel. The properties are listed in Table 3.2 for different temperatures. 


Temp. 

CC) 

Thermal Conduc. 
(W/m^C) 

Specific Heat 
(J/Kg"C) 

Coeff. of thermal exp. 
(lO-^/^C) 

Young’s Mod. 
(GPa) 

Yield Str. 
(MPa) 


37.8 

449 

11.6 

212 

425 


38.3 

468 

11.8 

208 

415 

100 

38.8 

488 

12.1 

205 

— 

200 

38.6 

526 

13.3 

200 

— 

mM 

37.5 

566 

13.8 

195 

ESHH 


Table 3.2 Material Properties for SA 508-Class 3 steel. 


3.4 Solid Modeling and FE Analysis 

The general purpose Finite Element code, ANSYS 5.4, has been used to perform the stress analysis 
of the structure. At first, a nozzle-vessel junction with the given geometry has been modeled, 
without fillets, and stress analysis has been performed without the presence of crack, under a 
pressure loading and taking linear elastic material behaviour. The corresponding solid model 
generated in ANSYS is shown in Figure 3.2. 
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Since the structure is symmetrical about the plane containing the axes of both the cylinders 
(plane 1, as shown in Figure 3.2), and also about the plane containing the axis of the no 2 zle and 
perpendicular to the axis of the vessel (plane 3), it is sufficient to model a quarter section of the 
no 2 zle. Also, for simplicity, only a quarter section of the vessel has been modeled. This is valid 
under the assumption, that any radial plane of a cylinder far away from local geometrical 
discontinuities will have displacement only in the radial direction. 



Figure 3.2 Solid model of nozzle-vessel intersection. 


To decide upon the optimum lengths of the nozzle and the vessel that should be taken for 
analysis, a number of different models with different nozzle and vessel lengflis are analyzed, and, 
hoop and longitudinal stresses are plotted along the lengths for both the cylinders. The criterion for 
the optimum lengths is that, the above stress plots should reach the corresponding values as 
obtained from the Lame’s thick cylinder formulae. This is to ensure that the longitudinal 
boundaries are sufficiently far away from the junction, so that, they are not affected by the local 
discontinuity stresses, and uniform longitudinal tractions can be applied on them as boundary 
conditions. 
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Based on the above, it has been found that 1.5 m vessel length (on one side of the junction) 
and 0.5 m nozzle length are sufficient to meet the above criterion. The corresponding FE model is 
shown in Figure 3.3. 



Figure 33 FE Discretization of the model. 


The model has been discretized using a combination of 20-noded brick elements and 10- 
noded tetrahedral elements, with 13-noded pyramid transition elements in between. The model 
consists of 807 hexahedral elements and 2337 tetrahedral and transition elements, i.e., a total of 
3144 elements, and 8251 nodes. The value of Young’s modulus (£) used is 210 GPa, and that of 
Poisson’s ratio (l>) is 0.3. 

3.4.1 Load and Boundary Conditions 

The load and boundary conditions applied to the discretized model (with reference to Figure 3.2) 
are as follows: 
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a) With reference to Figure 3.2, Planes 1 and 3 are restricted to move in the respective 
normal directions, since both of them are planes of symmetry. 

b) Uniform internal pressure of magnitude 15 MPa is applied on the inner surfaces of both 
the cylinders. 

c) Plane 2 is restricted to move in the normal direction, since it is a radial plane to a 
cylinder loaded by internal pressure, and therefore cannot move tangentially. 

d) Longitudinal stresses are applied at the ends of both the cylinders, as calculated from 
the following relation: 


C7 = 



(3.1) 


where, R, and represent the inner and outer radii of the cylinders, and P is the 

internal pressure. This boundary condition is based on the assumption that the ends of 
both the cylinders are effectively capped by rigid closures. 


3.5 Results and Discussion 

The results of stress analysis are presented in the form of plots in Figures 3.4-3.7. In these figures, 
hoop and longitudinal stresses for the cylinders as calculated on the plane 1 are plotted along the 
lengths of the respective cylinders. Figure 3.4 shows the variation of stresses on plane 1 in the 
hoop and longitudinal directions to the vessel, along the edge which corresponds to the inner radius 
of the vessel (Path A-B). Figure 3.5 shows the variation of the stresses along the edge which 
corresponds to the outer radius of the vessel (Path C-D). Similarly, Figures 3.6 and 3.7 show the 
variation of the stresses on plane 1 in the hoop and longitudinal directions to the nozzle, along the 
edges corresponding to the inner (Path A-E) and outer (Path C-F) radii of the nozzle respectively. 

As clearly seen from Figures 3.4-3.7, all the stresses become constant towafds the end, as 
expected. The values to which these stresses approach are compared to the corresponding values 
calculated from the Lame’s thick cylinder formulae, and are found to match well within 2%. Thus, 
the validity of the boundary conditions is established. 
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To gain further confidence in the results of the stress analysis, the nozzle-vessel geometry 
analyzed by Moini and Mitchell [19], has been analyzed on the same lines as used in the analysis 
of the geometry of interest in the present work, and the results are compared with those of Moini 
and Mitchell [19], The dimensions of geometry analyzed, load applied and properties used by 
Moini and Mitchell [19] are: R,= 5 in, i?^= 5.5 in, r^= 2.5 in, r/= 2.75 in, r,= r„= 0.25 in. 

Internal Pressure P=100 psi, E= 30x10* psi, u= 0.3. 

The results of the present analysis are found to match satisfactorily well with the results of 
Moini and Mitchell [19], as shown in Figures 3.8-3.15. 



Figure 3.4 Stresses on the inside of the vessel (Path A-B). 
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Figure 3.5 Stresses on the outside of the vessel (Path C-D). 
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Figure 3.6 Stresses on the inside of the nozzle (Path A-E). 





Figure 3.7 Stresses on the outside of the nozzle (Path C-F). 



Figure 3.8 Hoop stress on the inside of the vessel (Path A-B). 
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Figure 3.9 Longitudinal stress on the inside of the vessel (Path A-B). 
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Figure 3.10 Hoop stress on the outside of the vessel (Path C-D). 
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stress (psi) ► Stress (psi) 



Figure 3.14 Hoop stress on the outside of the nozzle (Path C-F). 
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Figure 3.15 Longitudinal stress on the outside of the nozzle (Path C-F). 


3.6 Closure 

The results of stress analysis of the nozzle-vessel intersection without crack have been presented. 
By the comparison of the present analysis with the work of Moini and Mitchell [19], it is 
concluded that, the partial modeling of the nozzle-vessel geometry in a judicious manner gives 
reasonably satisfactory results while keeping the computational effort minimum. 
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Chapter 4 


EFFECT OF NORMAL OPERATIONAL LOAD ON 
NOZZLE-CORNER CRACKS 


4.1 Introduction 

As introduced in section 1.1, a high stress concentration exists at the crotch region of a nozzle- 
vessel intersection. Moreover, because of the welded joint at the intersection, there always exists a 
possibility of initial flaw. Due to the high stresses, an existing flaw may grow up in to a surface 
crack through fatigue. Such cracks are called nozzle-comer cracks. Though the actual shapes of 
such cracks are not simple, it has been observed that nozzle-comer cracks grown by fatigue are 
usually almost quarter-circular or quarter-elliptical in shape. Hence, for the purpose of numerical 
analysis, such cracks can be approximated by quarter-circular or quarter-elliptical profiles. In this 
chapter, the effect of the normal operational load (constant internal pressure) on nozzle-comer 
cracks has been studied. 


4.2 Crack Configuration 

Since the present work is a safety assessment, the most crucial crack orientation has been studied. 
The face of the comer crack has been assumed to be located in the plane which contains the axes 
of both the cylinders (vessel and nozzle), so that, the hoop stresses of the vessel and the nozzle 
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cause the crack to be loaded in mode I (crack opening mode, i.e. load normal to the crack 
surface). Figure 4.1 defines the dimensional parameters used to characterize the crack. 

With reference to the Figure, the centre, O, of the quarter ellipse, is characterized by r, , the 
inner radius of the nozzle, and i?, , the inner radius of the vessel. The crack depth a, along the 
nozzle wall, characterizes the semi-major axis of the quarter-ellipse, and similarly, the crack width 
b , along the vessel wall, marks the semi-minor axis of the quarter-ellipse. The same definition 
holds good when the inner fillet is also modeled, in which case, point o does not physically exist. 



Vessel centre Sne 



Ri 


Figure 4.1 Crack Configuration and dimensional parameters. 


It may be worth mentioning here, that the dimension of the comer crack in the nozzle 
direction is naturally more than that in the vessel direction, because in this case, the SIF along the 
crack front is more uniform, and the natural tendency of the crack is to grow in such a manner that 
the variation of SIF along the crack front is least. 

The aspect ratio of the crack is defined as the ratio of the crack width (h) to the crack 
depth (a ), so that, it is always less than unity. 
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The location of any point P, on the crack front is characterized by the angle ^ as shown in 
Figure 4.1. 


4.3 FE Modeling and Analysis 

Again, ANSYS 5.4 has been used for the modeling and analysis of the nozzle-vessel intersection 
with comer crack. Due to the presence of the crack at the above described location, the S 5 nnmetry 
condition about the plane containing the axis of the nozzle and perpendicular to the axis of the 
vessel (Plane 3 in Figure 3.2), which prevails in the absence of the crack, is not prevalent now. 
Therefore, in this case, a semi-section of the nozzle has been modeled along with a quarter section 
of the vessel, rather than a quarter section of the nozzle, which is the case in the absence of the 
crack. A vessel length of 3 m and a nozzle length of 0.5 m have been taken in the model (Figure 
4.2). 

Since the crack is in the crotch region, it is realized that the inner fillet may have a 
significant effect on the crack. But the outer fillet is still far away from the crack and is unlikely to 
have any significant effect on the stress state near die crack. That is, the inner fillet falls within the 
region of interest, and hence, it has been included in the model. 



Figure 4J, SoUd model of nozzle-vessel intersection with comer crack. 
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The details of the modeling of the crack are shown in Figure 4.3. The dimensions of the 
crack are; Crack Depth = 80 mm, Aspect Ratio = 0.8. A small region above, and below the crack 
front (with a thickness of 5 mm in the direction normal to the crack front); have been modeled as 
separate regions, with the intention to mesh them with hexahedral elements (rather than tetrahedral 
elements). This is because, as it will be seen in Section 4.7, a structured mesh (not possible with 
tetrahedral elements) is required near the crack front, for the purpose of evaluation of SIF. Again, 
each of these smalt regions has been further divided in to five sub-regions by partitions orthogonal 
to the crack front. This is done to ensure the meshing of these volumes in such a way that the 
element edges in the plane of the crack face (the de-bonded surface above the crack front is 
referred to as crack face, as shown in Figure 4.3) are orthogonal to the crack front. This is required 
(as will be seen in Section 4.7) for a better accuracy in the evaluation of SIF. 


Vessel Side 



Nozzle Side 


Figure 4.3 Enlarged view of the near-crack region. 
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The above model (Figure 4.2) is then meshed with a combination of 20-noded hexahedral 
elements and 10-noded tetrahedral elements. The region near the crack front has been meshed with 
160 hexahedral elements, with 20 elements along the crack front, 2 elements above and 2 elements 
below the crack front, and 2 elements in the direction perpendicular to the crack face (as shown in 
Figure 4.4). The total number of elements in the discretized model is 15386, comprising of 2348 
hexahedral elements and 13038 tetrahedral and transition elements. The total number of nodes is 
32804, which amounts to atotal of 98412 degrees of freedom. 

There are two main reasons for this manifold increase in the number of elements required 
for the discretization, as compared to the case of stress analysis without crack. Firstly, a relatively 
finer mesh has been used in the crotch region, to effectively capture the high stress gradients 
existing near the crack. Secondly, the additional portion that has to be modeled because of the non- 
symmetry about Plane 3, due to the presence of the crack. 



Vessel Side 


Nozzle Side 


Figure 4.4 FE discretization of the near-crack region. 
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A fine mesh near the crack is essential, however, since the effect of a crack is local, it is 
possible that the symmetry condition about Plane 3 is actually not disturbed much, and the 
symmetry boundary condition is still applicable to the plane. 

To test the validity of the above realization, stress analysis has been performed with the 
full model (Figure 4.2), and also with a half model (as shown in Figure 4.5), assuming symmetry 
about Plane 3, and the results are compared. The total number of elements in the half model is 
9971, comprising of 1254 hexahedral elements and 8717 tetrahedral and transition elements, and 
the total number of nodes is 20484. The material properties used in the analyses are; Young’s 
modulus, JE = 210 GPa, and Poisson’s ratio, d = 0.3. 



Plane 2 


\ 
Plane 1-^ 


Plane 3 

I — * 


Figure 4.5 Half model with definitions of paths and coordinate axes. 
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4.3.1 Load and Boundary Conditions 


The load and boundary conditions applied during the stress analyses, with reference to 
Figure 4.5, are as follows: 

Full model: 

1 ) Symmetry boundary condition is applied to the Plane 1 , excluding the crack face. 

2) Plane 2 is restricted to move in the normal direction. 

3) Internal pressure of magnitude 15 MPa is applied on the inner surfaces of both no 22 ;le and 
vessel, and also on the crack face. 

4) One of the longitudinal ends of the vessel is fixed in the normal direction and longitudinal 
traction is applied on the other, as described in section 3.3.1. 

5) Traction boundary condition is applied to the longitudinal end of the nozzle, as per section 
3.3.1. 

Half model: 

1 ) Symmetry boundary condition is applied to the Plane 1, excluding the crack face. 

2) Plane 2 is restricted to move in the normal direction. 

3) Symmetry boundary condition is applied to the Plane 3. 

4) Internal pressure of magnitude 15 MPa is applied on the inner surfaces of both nozzle and 
vessel, and also on the crack face. 

5) Traction boundary conditions are applied to the longitudinal ends of the vessel and the 
nozzle, as per section 3.3.1. 

4.3.2 Comparison of Results for Full and Half Models 

To compare the results of stress analyses for full and half models, two paths have been chosen on 
Plane 3, corresponding to the inner and outer radii of the vessel. Path G-H and Path I-J (as shown 
in Figure 4.5), and various stresses are plotted along these paths, for both full and half models and 
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are compared in Figures 4.6-4,l 1. In Figure 4.12, the variation of the crack opening stress (stress in 
X-direction) along the crack front has been compared for full and half models. 

As it is obvious from Figures 4.6-4.11, the state of stress on Plane 3 is more or less the 
same for full and half models. This indicates that Plane 3 does not really feel the presence of the 
crack much, and therefore it can still be considered as a plane of symmetry. 

It is observed from Figure 4.12, that the crack opening stress, as estimated by the half 
model, is slightly on the higher side consistently, the maximum difference being less than 2%. This 
is expected, since, the assumption of symmetry about Plane 3 is equivalent to the presence of 
another comer crack located symmetrically on the other side of Plane 3. 

It is concluded that, by the assumption of symmetry about Plane 3, a very little error creeps 
in to the solution, and that too, on the conservative side, and hence the half model can be used for 
stress analysis of the nozzle-vessel intersection with comer crack. 

Finally, quarter-point elements (as discussed in Section 4.7) have been incorporated near 
the crack in the half model, to simulate the square root singularity that exists near a sharp crack 
(Section 4.5), and stress analysis has been performed. 20-noded brick elements with the quarter- 
point configuration of the mid-side nodes (adjacent to the nodes S), as shown in Figure 4.23(c), 
have been used as singular elements. The results of stress analysis of the half model after 
incorporating the singular elements are discussed in Section 4.4. 


4.4 Stress Analysis Results and Discussion 

The results of the stress analysis after incorporation of singular elements are presented in the form 
of plots in Figures 4.13-4.18. These plots show the variation of stresses in different directions 
along the Paths A-B and A-E (defined in Figure 4.5). It is clearly seen that the stresses at the crack 
tip become very high as compared to the stresses in the surrounding. Both the paths fall along the 
element edges of the singular elements near the crack front, and hence, square root nature is 
expected in the singularity. Figures 4.19 and 4.20, which show closer views of the stresses near the 
crack tips, demonstrate the square root nature clearly. Square root singularity is well achieved 
within the first element (which is quarter-point element) beyond the crack front. 
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Figure 4.8 Stress in Z-direction along path G-H 
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Figure 4.9 Stress in X-direction along path 1-J 
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Figure 4.10 Stress in Y-direction along path I-J. 



500 1000 
Distance (mm) ► 


Figure 4.11 Stress in Zrdirection along path hJ. 
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4.5 Stress Intensity Factor 


U V, ni 231 that the stress field at any arbitrary point, Q, in the vicinity of the crack 
It can be shown [22, 23] that the the case of mode I loading, in an 

up T (Figure 4,21). for isouopic, linear elasuo matenal, for fte case of m 

infinite flat plate »ith a through-thickness crack of length 2 n, is grvenby. 
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(4.1c) 
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where, <t is the far field stress; r and 6 are the polar coordinates of the point Q, with the origin 
stationed at the crack tip T, and 9 being measured anticlockwise from the extended crack line. 



Figure 4.21 Through-thickness Crack in an infinite flat plate. 



crack front and crack faces 


Figure 4.22 Stress field around an embedded elliptical flaw. 
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The coordinate axes, AT, and Xj , being along, and perpendicular to the extended crack line, 
respectively. Here, denotes the component of stress that acts in the direction of the j"’ 

coordinate axis on a plane whose normal is along the coordinate axis. 

For a thin plate (plane stress case), the other stress components are negligible. In case of a 
thick plate (plane strain case),a 33 = l>(<Th where, v is the Poisson’s ratio of the material; 

and the remaining two stress components are negligible. 

It can be seen from Equations 4.1, that the distance (r ) between the crack tip and the point 
occurs in the denominator under a square root sign, in all the stress components. If r becomes very 
small, the stress components rise very sharply and tend to become infinite as r — >■ 0 . This means 
that the stresses become singular at the crack tip, and the particular kind of singularity encountered 
here is known as square root singularity because of the xj-Jr nature of the stress field. 

Also, it can be observed from Equations 4.1, that the quantities cr and a always coexist as 
<7 in all the equations of the stress field. Apart from these quantities, the other variables 
(rand^) occurring in the field equations are merely geometrical variables defining the location of 
a point in the vicinity of the crack tip. This indicates that the product crV^ completely 
characterizes the state of stress existing in the vicinity of the crack. Hence, it is advantageous to 
recognize the product as a single parameter which can characterize the state of the crack 
completely. 

The idea was first proposed by Irwin in 1957, who defined the product as a new parameter, 
stress intensity factor, which is formally defined as: 

K, = (Itct) cr^i (r , ^ = 0) , as r 0 , (4.2) 

where, the subscript / denotes mode I loading. It can be verified that by substituting Equation 
4.1b in the formal definition (Equation 4.2), one can obtain K, = cr4m . 

The field equations can now be conveniently re-written in terms of the stress intensity 
factor K, . The corresponding displacement field in the vicinity of the crack tip, in terms of K, , is 
given by: 
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«3=0, (4.3c) 

where, E is the Young’s modulus of the material; and Ar = 3-4L> for plane strmn, and 
(3 — o)/ (1 + v) for plane stress. Here, m, denotes the displacement along the i'* coordinate axis. 


4.6 SIF for Three-Dimensional Crack Problems 

Real life crack problems are usually three dimensional, and involve part-through cracks, rather 
than through-thickness cracks. Part-through cracks (embedded or surface cracks) usually have 
curved crack fronts, which are modeled as elliptical, semi-elliptical or quarter-elliptical in the 
literature of fracture mechanics. 

Kassir and Sih [24] have demonstrated that, for the case of an embedded elliptical flaw, the 
singular stresses in the vicinity of any arbitrary point on the flaw border, in a plane mutually 
perpendicular to the flaw border and the crack faces, take the same functional form as the Irwin 
field equations for the two-dimensional problem. 

Figure 4.22 describes the above statement pictorially. Consider a cracked surface bounded 
by an elliptical crack front as shown in Figure 4.22. To understand the state of stresses around the 
crack profile, consider any point P (crack tip) on the crack front. Assume a plane cutting the crack 
faces such that, it is normal to the crack front at point P, and also perpendicular to the crack faces. 
Then, according to the work of Kassir and Sih [24], the stress field around the point P in the above 
plane is of the same form as that of the two dimensional problem of through-thickness crack in an 
infinite flat plate, and the stress and displacement components at a point Q (Figure 4.22) on the 
plane, in the vicinity of the point P are given by Equations 4.1 and 4.3 respectively, with the 
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definitions of the coordinate axes and the polar coordinates of point Q, as shown in Figure 4.22. 
This defines the SIF at point P. Thus, for a three-dimensional crack, there exists a continuous 
distribution of SIF over the crack front, rather than a single value. 


4.7 Numerical Evaluation of SIF 

The procedure for numerical evaluation of SIF has been described in detail by Hellen and Dowling 
[3]. The method makes use of the basic displacement field relations given by Equations 4.3. To 
evaluate the SIF at a given crack tip, the values of displacement (from the results of stress analysis) 
at nodes around the crack tip node are conveniently used in these equations, although first it is 
necessary to transform the displacement components from the global to the local directions as per 
Figure 4.22. 

Since the displacement field equations are valid in a plane mutually perpendicular to the 
crack front and the crack faces, it is necessary to have nodes in such planes corresponding to all the 
crack tip nodes. This is possible only if structured mesh is used near the crack front. With a 
structured mesh using quadratic brick elements, nodes exist along the 6 -=71 jl line, and also along 
the 6 = k line (if it is ensured that the meshing is such that the element edges in the plane of the 
crack face are orthogonal to the crack front). But it has been observed that the results of 9 
line are more accurate than those of ^ = Tijl line [3]. 

Denoting the component of displacement in the direction normal to the crack face by v . 
one can obtain, from Equation 4.3 b: 
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for 9 = n line. 
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Since two nodes lie along either line in the crack tip elements, two independent results for Kj are 
obtained using either of the Equations 4.4. However, Equations 4.3 do not include the 
homogeneous strain terms which occur in reality, such as a large hoop stress or thermal strains. 
Since the finite element model includes such effect as a linear function of r, a linear correction is 
required, to give Equation 4.4 as; 


v = C 


K,(r)y-‘ ^r 


Lr 


(4.5) 


where, C and L are constants to be eliminated. Using the subscript m for values at the mid side 
node and v for values at the comer node, and noting that = 4r„ (for quarter-point crack tip 
element), 


= 




(4.6a) 


and 


v„ =C 


K^(ry^+Lr, 


(4.6b) 


However, Equation 4,3b may still be used if an equivalent value K] is used such that: 

v = CKXy. (4.7) 

By inspection of Equations 4.6 and 4.7, use of such a K'j requires die following conditions to be 
true: 


and 



(4.8a) 


(4.8b) 
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By suitiiblc manipulation of Equations 4.8, and using , one obtains! 

(4.9) 

Equation 4.9 can be used to evaluate the SIF values at the crack tip nodes using the displacement 
field obtained from stress analysis. Improved results are obtained [3] by the use of Equation 4.9. 

Also, it is important to simulate the square root singularity existing at the crack front, 
during the stress analysis, to get reasonable results. One of the simplest and most popular methods 
of simulating singularity is the use of quarter-point elements near the crack. An element of this 
kind [25], shown in Figure 4.23(a), introduced almost simultaneously by Henshell [26] and 
Barsoum [27] for 2D problems, develops by a simple shift of the mid-side nodes (adjacent to the 
comer node at which singularity is desired, indicated by ‘S’ in the Figure) in quadratic, 
isoparametric quadrilateral elements, to ihe quarter point locations. 



(a) 2D quarter point (b) 2D side coHapsed 

eiement quarter point element 



(c) 3D quarter point (d) 3D face coHapsed 

element. quarter point element. 


Figure 4.23 Quarter point Eientents for simulation of singularity. 
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Singularity develops at the node S because, in the particular quarter-point configuration of the mid- 
side nodes, the deteiminant of the Jacobian becomes zero at the node S. It can be shown that along 

the element edges, the derivatives dufdx (or strains), and hence stresses, vary as l/ -Jr , where, r 
is the distance from the comer node S, at which the singularity develops. 

A better element for the purpose, as shown in Figure 4.23(b), is the side collapsed quarter- 
point element [24] proposed by Hibbitt [28], in which, square root singularity is achieved in any 
radial direction emanating from the comer node S. 

Although singularity is not well modeled on lines other than element edges in the 
quadrilateral quarter-point element, as it is done in case of the side collapsed element, yet, 
reasonably good results are achievable [25] with such elements. 

All the above ideas can be extended to 3D [29], and quarter-point elements can be 
developed from second order 3D elements, as shown in Figures 4.23(c) and (d), for simulation of 
singularity in 3D crack problems. 

The above described technique has been used to evaluate the SIF values for nozzle-comer 
cracks in the present work. SIF values have been evaluated along the crack profiles for two 
different crack configurations: 

1) Crack- 1 : Crack Depth = 80 mm. Aspect Ratio = 0.8; representing medium to large cracks. 

2) Crack- II : Crack Depth = 35 mm. Aspect Ratio = 1 .0; representing shallow cracks confined 

within the fillet region. 

Stress analysis has been performed for both the crack configurations using quarter-point elements 
near the crack front, with the same load and boundary conditions as described in Section 4.3.1. SIF 
values for Crack- 1 have been evaluated along the profile, for ^ = /r line. 

For Crack- II , it has not been possible to mesh the region near the crack with element edges 
in the crack plane orthogonal to the crack front throughout the profile, because of the intersection 
of the crack front with rounded edges. Therefore, for Crack- II , SIF values have been evaluated for 
B = 7tl2 line. 

The plane stress relations have been used at the end points of the profile (which lie at free 
surfaces), and plane strain relations have been used at all inner points, for both the cracks. 
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4.8 SIF Results and Discussion 


The variation of SIF along the crack prbflle for Cracks-I and II. as evaluated by the method 
described m Section 4.7, have been shown in Figures 4.24 and 4.25 respectively. In these figures, 
the SIF values at the crack tip nodes have been plotted against the angular position (angle ^ as 

shown in Figure 4.1) of the nodes. The angle 0® corresponds to the vessel-side end point of the 
crack front, while the angle 90 corresponds to the nozzle-side end point of the crack front. 

It is seen that, for Crack- 1 , the SIF has higher values near the free surfaces (except in very 
thin regions near the surfaces), and tend to decrease towards the inner side, with a (local) minimum 

value at around ^ — 45 . The average SIF over the profile is about 848 MPa-J mm , with a 
percentage variation of SIF of about 21 % over the profile (defined as the percent ratio of the 
difference between the maximum and minimum SIF values over the profile excluding the end 
points, to the average SIF). The higher values of SIF near the end points may be expected, because 
these regions exist near the inner radii of the cylinders, where the hoop stress has relatively higher 
values. It should also be noted that the SIF is higher near the nozzle side as compared to that near 
the vessel side. This indicates that such a crack will become critical first on the nozzle end and 
later on the vessel end (for pressure loading only), and hence there is a chance for leak-before- 
break. 

It can be seen that the SIF falls off sharply in very thin layers near the free surfaces. This is 
not surprising, since very near to the free surfaces, plane stress condition exists. This tendency has 
also been verified from the work of Mohamed and Schroeder [9], wherein, it is clearly mentioned 
that the SIF actually drops near the free surfaces. 

In the case of Crack- II , which is completely contained within the fillet region, it is seen 
that the pattern of variation of SIF along the profile is entirely different. The SIF values are higher 
inside and drop very sharply near the free surfaces. It has been verified from the work of Heliot et 
al. [30], that the shape of the SIF curve for shallow cracks confined within the fillet region is 
actually similar to the one shown in Figure 4.25. However, the SIF is still higher on the nozzle side 

than on the vessel side. The average SIF over the profile is about 482 MPa-J mm , with SIF 
variation of about 25 % over the profile (excluding the end points). 
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4.9 Validation of SIF Results 


To test the validity of the SIF results obtained in the present work, and to get an estimate of the 
order of accuracy of the results, two test problems (one 2D and the other 3D), whose handbook 
solutions are available, have been solved. 

Also, the SIF results of the nozzle-comer crack problem as obtained from the present work 
are compared with the results obtained from an approximate empirical formula available in the 
literature. 


4.9.1 2D Test Problem 

A single-edge-notch-bend (SENS) specimen, as shown in Figure 4.26, has been analyzed using the 
8-noded quadrilateral elements, with quarter-point crack tip elements. The plane strain condition 
has been assumed. With reference to Figure 4.26, the dimensions of the specimen analyzed are: 

Span, S=lm, Width, W= 0.25 m, Thickness, B=lm (plane strain). Crack length, a/W = 0.5. 

The load applied, ? = 100 iV . Only half the model is analyzed, exploiting the symmetry of 
the problem. SIF is evaluated using the method described in Section 4.7. The value of SIF obtained 

is 2193.6 PaVm . 

The SIF is also calculated from the Handbook Solution for SENB specimen [22], which is 
as follows: 


K, = , where, /(0.5) = 2.663 . (4. 1 0) 

The value of SIF obtained from this relation is 2130.4 Pa^ . The error in the present result 
comes out to be + 2.97 % (conservative). 
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4.9.2 3D Test Problem 


This problem is taken from Murakami’s Stress Intensity Factors Handbook [31]. The problem is 
that of a quarter-elliptical corner crack emanating from a central hole in a finite flat plate subjected 
to uniform tension, as shown in Figure 4.27. The crack is located in the diametral plane of the hole, 
parallel to the loaded face of the plate. The dimensions and load used in the analysis are: 


Height of the plate, 2H = 2400mm , Width of the plate, 2W = 1400 mm , Thickness, r = 1 00 mm , 
Radius of the hole, /? = 1 00 mm , Crack Depth, a = 25 mm. Crack Width, b = 50mm. 

The load applied is a uniform mode I stress, a = 15 MPa as shown in the Figure. 

Only one-fourth of the geometry has been modeled (assuming quarter-symmetry) and 
analyzed, using 20-noded brick elements and incorporating quarter-point crack tip elements. SIF 
has been evaluated along the crack profile using the plane strain relations throughout. SIF is also 
calculated from the Handbook Solution [31], and the variation of both the SIF values along the 
crack profile have been compared in Figure 4.28. 

The Handbook Solution is reported to match with FEM results within ±5%. It shows a 
maximum difference of about ± 7% from the present FEM results, except for the first point. The 
difference in average SIFs is only 0.4%. From this comparison, it seems that the present FEM 
results are within ± 2% of the standard FEM results for the problem considered. 


4.9.3 Comparison with Approximate Empiricai Formula 

% 

The present results for a quarter-circular nozzle-comer crack are compared with the results 
obtained from the approximate empirical formula proposed by Mohamed and Schroeder [9]. The 
formula relates the SIF at a point on the crack front to the hoop stress existing at that point in the 
absence of the crack, employing the solution of an elliptical crack embedded in an infinite solid 
exposed to uniform tension, with some modifications. This is based on the assumption that the 
uncracked stress distribution reflects the influence of all the geometrical variables affecting the 
Stress Intensity Factor. 


57 



The SIF at any point P, at an angular position <f> , on the front of a noTzle-comer crack with 
depth = and width = 6^ , is given [9] as: 




G, 


■yjm 


E(k) 


.2 . . 


sin‘‘^+-|-cos ^ 
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(4.11) 


where, a is the crack length at ^ = 45", is the magnification factor for the front free surfaces 
which is given as a function of (j) in the form of a curve for different values of aspect ratio, cr^ is 
the hoop stress at point P in the absence of the crack, and, E{k) is the complete elliptical integral 
of the second kind: 
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From Equations 4. 1 1 and 4. 12, the SIF for a quarter-circular crack is given by: 





(4.12) 


(4.13) 


Figure 4.29 shows the variation of SIF along the crack front of an 80 mm deep quarter- 
circular nozzle-comer crack as obtained from the present work, and also from the empirical 
formula (Equation 4.13). In the present SIF results, plane stress relations have been used at the end 
points of the crack front, and plane strain relations have been used at all interior points. 

The above comparison reveals that the results of the present work are matching, with a 
maximum difference of about 20 %, with the results of the empirical formula. The difference 
between the average SIFs is about 5 %. 

It is reported by Mohamed and Schroeder [9], that their results match with the results of 
FEM with a maximum deviation of the order of 16 %. In lieu of the above, it seems that the present 
results match within 4 % with the standard FEM results for nozzle-comer cracks. 
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Figure 4.28 Comparison of SIF for 3D test problem. 



Figure 4.29 Comparison of SIF for 80 mm deep qtr-circular nozzle-corner crack 
with the results obtained from empirical formula [9|. 
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4.10 Parametric Study 


To understand the effect of different geometrical and crack parameters on the SIF of nozzle-comer 
cracks, parametric study has been performed by varying the parameters in the neighbourhood of 
their given values. Figures 4.30-4.33 show the variation of SIF along the profile of nozzle-comer 
cracks, for different values of crack and geometrical parameters. Since the analysis is linear with 
respect to the load (internal pressure), the load has not been taken as a parameter in the parametric 
study. All the results presented here correspond to an internal pressure of 15 MPa. 


4.10.1 Normalization of the Parameters 

The different geometrical parameters involved in the problem are: R^ (inner radius of the vessel), 
R 2 (outer radius of the vessel), r, (inner radius of the nozzle), (outer radius of the nozzle), and r 
(inner fillet radius). 

For the purpose of parametric study, these parameters have been taken as: D (mean 
diameter of the vessel), T (thickness of the vessel), d (mean diameter of the nozzle), t (thickness of 
the nozzle), and to reduce the number of independent variables, the inner fillet radius (r) has been 
taken equal to the same fraction of the thickness of the vessel wall, as it is ih the given geometry. 
This is in accordance with the thumb rule of design of pressure vessels [32] to take the inner fillet 
radius as some fraction of tlie thickness of the vessel wall. 

Now, these four geometrical parameters have been combined into three dimensionless 
groups in the following manner, there by reducing the number of independent geometrical 
variables to three; 


£ L R. 

D’ T' T 

The above choice of dimensionless groups has been guided by the work of Mohamed and 
Schroeder [9], wherein, approximate empirical formula for SIF of nozzle-comer cracks has been 
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given m terms of the groups ^ and — . The group ~ has been added to the set for the sake of 

completeness, since, four independent variables can at most be reduced to three independent 
dimensionless groups. 

The crack parameters are: crack depth {a) and aspect ratio (b/a). The crack depth is 

normalized with respect to the thickness of the vessel wall as- — 

■ T' 


The SIF has been normalized with respect to the internal pressure (P), and the average 


crack length, a 


a+b 


,as: 




4.10.2 Discussion 

In Figures 4.34-4.40, variation of SIF has been plotted with one parameter at a time, keeping all 
others constant. Figure 4.34 shows the effect of crack aspect ratio on the variation of SIF over the 
crack profile. Three important conclusions can be clearly drawn from the Figure. Firstly, the 
average SIF considerably reduces, as the crack becomes more and more elliptical, i.e. as the aspect 
ratio is reduced from unity. Secondly, the percentage variation of SIF from the mean value reduces 
considerably as the aspect ratio decreases. That is, the variation of SIF becomes more and more 
uniform. Thirdly, the point of minimum SIF shifts towards the nozzle side as the aspect ratio 
decreases. 

Figure 4.35 shows the variation of normalized SIF over crack profiles for different aspect 
ratios. It is seen that the normalized SIF plots for different aspect ratios are closer to each other as 
compared to the SIF plots for different aspect ratios (Figure 4.34). Interestingly and fortunately, it 
so happens that the normalized average SIF remains almost constant with respect to aspect ratio, as 
shown in Figure 4.36. This means that the variation of normalized average SIF with respect to any 
other parameter will not depend on the value of aspect ratio. 

Normalized average SIF has been plotted against the normalized crack depth in Figure 4.37 
for a given Geometry. The decreasing nature of the curve has been verified form the literature [9]. 
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Figure 4.38 shows the variation of normalized average SIF with respect to d/D. The 
increasing trend of the curve is expected, since, a larger d/D ratio signifies a more severe 
discontinuity in the vessel. 

Variation of normalized average SIF has been plotted against the ratio t/T in Figure 4.39. 
Figure 4.40 shows the variation of normalized average SIF with respect to the ratio D/T. The 
increase in normalized average SIF with increase in D/T is expected, due to the fact that the hoop 
stresses in the vessel increase with the increase in D/T (this can be verified from Lame’s thick 
cylinder formula). The almost linear nature of the curve may also be attributed to the almost linear 
dependency of the hoop stresses on the ratio D/T (this can also be shown using the Lame’s thick 
cylinder formula). 



Figure 4.30 Variation of SIF along the profiles of cracks with different depths for the given Geometry. 
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Figure 431 Variation of SIF along the profiles of cracks with different depths 
for three different Geometries, 



Figure 4.32 Variation of SIF along the profile of an 80 mm deep crack for two different Geometries. 
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Figure 4.33 V&riation of SIF along the profile of an 80 mm deep crack for three different Geometries. 



Figure 4.34 Variation of SIF along the profiles of cracks with a given depth for three different aspect ratios. 
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Figure 4.36 Variation of normalized average SIF with respect to crack aspect ratio for a given 
Geometry and crack depth. 
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Figure 4.37 Variation of normalized average SIF with respect to normalized crack depth 
for a given Geometry. 



Figure 4.38 Variation of normalized average SIF with respect to d/D. 
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4.11 Closure 


Stress analysis of the no22le-vessel intersection with comer crack has been performed using 
singular elements near the crack. Stress Intensity Factors have been evaluated along the crack 
profile and the methodology adopted has been verified. Finally, parametric study has been 
performed to understand the effect of different parameters on the SIF. 

The SIF of nozzle-comer cracks has been found to bear simple relationships with most of 
the parameters individually, indicating the possibility of development of more accurate empirical 
relations for its calculation, as compared to those already available in the literature. 
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Chapter 5 


EFFECT OF LOCA ON NOZZLE-CORNER CRACKS 


5.1 introduction 

As already introduced in section 1.1, LOCA refers to a loss of coolant accident in a nuclear reactor 
plant. Fracture failure of reactor pressure vessel may occur during a LOCA if the prevailing 
circumstances are adverse. Three conditions are necessary for such an incident to occur: 

1) A large upward shift in the nil ductility transition (NOT) temperature due to a combination 
of nuclear irradiation during service and the presence of high copper and nickel content in 
the vessel welds. 

2) The existence of an initial flaw on the inner surface of the vessel. 

3) A severe over-cooling transient caused by the splashing of cold water on the inner surface 
of the vessel by the emergency core cooling system during a LOCA. 

The over-cooling transient involves a sharp change in the temperature of the inner surface of the 
vessel. Simultaneously, the internal pressure would be expected to fall from the operating pressure 
to the atmospheric pressure. This typical thermal load is termed as “Thermal Shock” in the nuclear 
industry. If the pressure remains high during the transient, or falls off slowly, the event is termed as 

“Pressurized Thermal Shock” (PTS). 
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Sudden cooling of the inner surface of the vessel causes sudden contraction of the inner 
layers, but the outer layers contract slowly. This produces high tensile stresses in the hoop 
direction on the inner side of the vessel (and nozzle). Hence, a flaw existing at the inner side of the 
vessel (or nozzle-vessel comer) with an orientation which causes hoop stress to load it in mode I, 
can be severely affected. 

Albeit, LOCAs have occurred during reactor operation, for instance, the Three Mile Island 
II, because all three of the conditions mentioned above were not simultaneously satisfied, no 
catastrophic fractures of nuclear plant pressure vessels have been experienced. 

Nonetheless, the problem is clearly significant and must be given special attention. 


5.2 Problem Definition 

The present intention is to study the effect of a given PTS transient on a quarter-elliptical comer 
crack existing at the crotch region of the given nozzle-vessel intersection. Since the temperature of 
the inner surface of the vessel and the nozzle will vary with time during the transient, the 
temperature distribution in the wall of the nozzle-vessel intersection will be a function of time. 
Moreover, internal pressure is also time-varying. This causes time-varying thermal and mechanical 
stresses and, in turn, time-varying stress intensity factors. 

It is intended to study the variation with time, of SIF along the crack front during the given 
PTS transient, and compare the SIF values with the material fracture toughness for crack initiation 
to estimate whether the crack becomes critical or not. Also, it is intended to identify the most 
critical time, and location on the crack front. 

The geometry of the nozzle-vessel intersection considered here is the same as specified in 
section 3.1 . The crack considered is characterized by: 

Crack Depth, ajT = 0.5 
Aspect Ratio, i/a = 0.8 

The PTS transient is given as a thermal and mechanical load history on the inner surface of 
the vessel and the nozzle, as tabulated in Table 5.1 and also shown graphically in Figure 5.1 . 
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Time (s) 

Temperature (“C) 

Pressure (bar) 

0 

310 

170 

■■ 

230 

150 

1500 

60 

6 

8000 

60 

6 


Table 5.1 Given PTS transient 


The table (and the figure) specifies the value of temperature at the inner surface of the vessel and 
the no 2 zle, and the magnitude of the internal pressure at different time points. The intermediate 
values will be obtained by linear interpolation. The transient is of duration of 1500 s (25 minutes) 
only, but the load history is given up to 8000 s (about 2 Hrs. and 15 minutes). This is done to reach 
as close to the steady state as possible, after the PTS. 



Figure 5.1 Given PTS transient. 
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5.3 FE Modeling and Analysis 


ANSYS 5.4 has been used to perform the transient thermal -stress analysis of the nozzle-vessel 
intersection with comer crack. The geometry has been modeled along with crack and is meshed 
with 18729 elements comprising of 1254 hexahedral elements, and 17475 tetrahedral and 
transition elements. The total number of nodes in the mesh is 32318. Two layers of hexahedral 
elements have been used on either side of the crack front with 20 elements along the crack front. 

First, a static heat transfer analysis is done to obtain the temperature distribution existing at 
time, t = 0, which represents the normal operating condition of the vessel. Then, a transient heat 
transfer analysis has been performed over the entire time domain of 8000 s, to obtain the 
temperature distribution at a number of time points. 

The given transient temperatures are specified on the inner surface of the vessel and the 
nozzle. Convective boundary condition is specified on the outer surface of both the cylinders. The 
convective film heat transfer coefficient at the outer surface of a cylinder has been taken from the 
Heat and Mass Transfer Data Book by Kothandaraman and Subramanyan [33] and is given as; 

;i = 1.3l(Ar®“') CmWlm^K), (5-1) 

for a turbulent boundary layer, where, ^T is the difference between the temperature of the surface 
and the bulk temperature. Turbulent boundary layer is chosen because for a laminar boundary 
layer, the heat transfer coefficient also depends upon the dimensions of the cylinder, whereas it is 
intended to use some uniform reference value for both the cylinders. Equation 5.1 is valid for a 
vertical cylinder and there is a different relation prescribed for a horizontal cylinder, but the 
difference is not much. Also, h is not a very strong function of AT. In lieu of the above, h is 
calculated from Equation 5.1 for an assumed average value of Alequal to 150‘’C (an atmospheric 
temperature of 25 “C is assumed). This value (« 7 Wlrn^K), is used as the average heat transfer 
coefficient for the outer surfaces of the vessel and the nozzle. All the other surfaces of the structure 
are assumed to be insulated. 

A total of 15 time points mre chosen among those at which temperature distributions have 
been determined. The immediate layer of elements on either side of the crack front is made 
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singular by shifting the mid nodes adjacent to the crack tip nodes, to quarter point locations. Stress 
analysis has then been performed at each of the 15 time points chosen. 

1 he internal pressure calculated by linear interpolation at the respective time is applied on 
the inner surfaces of both the cylinders and also on the crack face. The displacement and traction 
boundary conditions used are similar to those described in section 3.3.1. 

The material properties used in the analyses are taken from Table 3.2 at 100 °C assuming 
that it represents an average temperature during the transient. The other properties used in the 
analyses are: Poisson’s ratio, v = 0.3, and density, p = 7833 Kg/m^. 

Finally, SIF values have been evaluated (using the plane strain relations throughout the 
crack profile) along the crack front, for each of the 15 time points for which stress analyses have 
been performed. 


5.4 Results and Discussion 

The variation of temperature along the crack front has been plotted for different time points in 
Figures 5. 2-5.4. As seen in Figure 5.2, the temperature along the crack font is almost constant at 
time t = 0 s. With the outbreak of the transient, the temperatures at the end points (angular 
positions of 0® and 90®) of the crack front (which correspond to the intersection of the crack front 
with the inner surfaces of the vessel and the nozzle respectively), drop sharply, as seen in the plots 
for time t = 23 s and t = 50 s in Figure 5.2, whereas, the response of the points away from the «Ki 
points is slower. Except for a few initial seconds, the temperatures at the end points are lower and 
increase towards the middle. At all times, the maximum temperature along the crack front occurs 

at around 50®. 

The variation of SIF along the crack front has been plotted for different time points in 
Figures 5.5-5.7, along with the material fracture toughness for crack initiation (Critical Stress 
Intensity Factor, AT/,). The Critical Stress Intensity Factor for crack initiation for the material of 

the structure has been taken from ASME, Section XI, Appendix A [18], which gives a lower 
bound test data for K,^ as a function of temperature, as: 
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K,^ =33 2 + 2.81ex/;/0 02('r-/?7’j,,j, +100;; (in fcszV^), 


(5.2) 


subject to a maximum of 200 , where, T is the temperature (in"F), and is the 

reference NDT temperature {m‘’F). 

In the present analysis, RT^,^, has been taken equal to 0°C. Equation 5.2 attains the 
maximum value of 200 ksi-Jin at about 55 "C for this value of RT^cr , whereas the minimum 
temperature in the present analysis is 60 “C. Hence, the value of has been taken equal to 

200 ksi-Jln ( « 220 MPa4m « 6950 MPa^mm ) throughout the present analysis. 

Figure 5.5 shows that the SIF along the crack front is increasing with time, as expected. 
Initially, the SIF is higher towards the nozzle side (angular position = 90“), but as the PTS 
outbreaks, the SIF towards the vessel side (angular position = 0“) starts rising at a rate faster than 
that of the nozzle side. The SIF towards the middle of the crack front take still longer time to build 
up. The vessel side SIF soon becomes equal to the nozzle side SIF (at around t = 50 s) and 
surpasses it. It may be observed from Figure 5.6, that the point of minimum SIF is continuously 
shifting towards the nozzle side, as the time passes. 

The SIF keeps on increasing and finally reaches the peak at time, t = 1500 s (25 minutes), 
when it can be seen in Figure 5.7, that a small portion of the crack front (between angular positions 
of 0 “ and 6 “ ) near the vessel side has become critical (SIF value has exceeded Ki^). After t = 1500 

s, the SIF begins to fall off throughout the crack front. 

Figures 5.8-5. 10 show the variation of vessel side SIF, nozzle side SIF, and the average SIF 
over the crack front, with time, respectively. It may be observed from these Figures that the initial 
rate of rise of SIF is very high in all the three cases, but fall off continuously. The rate falls off 
more rapidly for the nozzle side SIF than for the vessel side SIF. After reaching the peak values, 
the SIF decay at an average rate much lower than the average rate of increase. The steady state is 

not reached even after 8000 s. 

It is apparent from Figure 5.8, that the crack front (towards the vessel side) is critical for a 
period of about 420 s (7 minutes), from time, t = 1340 s to t = 1760 s. The crack may or may not 
fail during this period, depending upon the degree of conservatism involved in the value o^K,^, 

and various other factors. 
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Figure 5.2 Variation of temperature along the crack front (Time points 1-5). 



Angular Position (degrees) — 


Figure 5,3 Variation of temperature along the crack front (Time points 6-10). 
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Figure 5.6 Variation of SIF along the crack front (Time points 6-10). 
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Figure 5.7 Variation of SIF along the crack front (Time points 11-15). 







Figure 5.9 Variation of average SIF with time. 


5.5 Closure 

Thermal-stress analysis of the nozzle-vessel intersection with a comer crack has been performed to 
study the effect of a LOCA on nozzle-comer cracks. A PTS appears to be much more cmcial than 
a normal pressure loading. The considered crack configuration seems to be unsafe, with a small 
portion of the crack becoming critical for a small duration of time during the PTS considered. The 
vessel side of the crack front has been identified to be the most critical location during a PTS, 
unlike a normal pressure loading, during which, the nozzle side has a higher value of SIF. This 
suggests that apparently there is no possibility of Leak-before-break during a PTS (as far as 
nozzle-comer cracks are concerned). However, the possibility of crack arrest needs to be assessed. 
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Chapter 6 


ESTIMATION OF STABLE CRACK GROWTH 


6.1 Introduction 

In the previous chapters, it has been attempted to estimate whether a crack is critical or not under a 
given set of loading conditions. When a crack becomes critical, die component becomes 
susceptible to catastrophic failure, but the crack growth takes place in a stable manner much before 
the crack becomes critical. Thus, to ensure a safe operation, it is important to estimate the amount 
of stable crack growth during a given load cycle, so that the safe life of the component may be 
accurately predicted, by estimating the number of load cycles required to grow the crack to a 

critical size. 

This chapter is an attempt towards the estimation of stable crack growth in a nozzle-comer 
crack using the cohesive zone technique. It is a numerical technique in which, the surface along 
which the growth of crack is expected, is modeled as a cohesive zone with a different matenal 
constitutive relation which is capable of simulating the softening, and ultimately failure, of the 
material. The formulation of the cohesive zone model is such that it can be conveniently 
implemented in a general purpose FE code through a user subroutine. The use of a cohesive zone 
model is inviting because the solution may be treated as a non-linear boundary value problem with 
a time dependent boundary without recourse to a crack growth criterion. 

The formulation of a 3D cohesive zone model, given by Foulk et al. [34] has been 
presented. The method is applied to a Double Cantilever Beam (DCB) specimen and a nozzle- 
comer crack and the results of crack growth analyses have been presented. 
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6.2 Constitution of the Cohesive Zone 


Cohesive zone models are employed to predict new surface areas evolving in a body, p as shown 
in Figure 6.1. By using constitutive models that are non-convex, it is possible to predict the 
creation of new surface area along any predetermined plane during the evolution of the problem. 

The constitution of a cohesive zone that is neither history dependent nor rate dependent is 
assumed to be of the following form [35]; 



where, the traction components, 7], and T, are coupled to both normal and tangential crack 
opening displacements, [m„], [m^J and [mJ respectively. Also, o-,^is a material property of the 

cohesive zone that is a measure of bond strength, and X is the normalized quantity coupling normal 
and tangential behavior; 


1 


kl' 


+ 


kl 








(6.4) 


where, and are material properties that are length scales associated with debonding, and 

<x ^ and cc^ are material properties relating shear to normal strength. It is furthermore assumed that 
the cohesive zone is fully debonded when ^ > 1 , so that in this case = — = 0 ^ It is to be 

noted that the displacement field is transformed to the local coordinate system shown in Fig 6.2. 
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Figure 6.1 Solid body with fixed and moving boundaries (cracks). 



Figure 6.2 Local coordinate system defined for each element. 
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T„ is depicted graphically in Fig 6.3 for the case of normal separation. In this case the 
maximum, , is reached when Ji= 1/3. If the cohesive zone is loaded beyond A= 1/3 and 
unloaded, the Tvergaard model will yield an increasing traction for a decreasing displacement. 
Because this is not physically realistic, the amended model as described by Tvergaard [35] is used. 
If /I > 1/3 and X<0, the following traction relationships are assumed to hold good: 






UJ 

^maxl c 


(l-2A + A') 



(6.5) 

( 6 - 6 ) 

(6.7) 


where is the maximum value of ^ before unloading. These relations remain valid 
until X > (reloading). 


6.3 Finite element Formulation 

The finite element formulation utilized herein is semidiscretized in that the displacement is 
modeled using shape functions that do not vary with time, whereas the time dependence m the 
problem is accounted for by using a finite difference (forward Euler scheme) approximation in 
time. Nonlinearity can be accounted for by a variety of iterative techniques. Herein, it is assumed 
that the Newton-Raphson method is used. 

6.3.1 Variational form of the boundary value problem 

The principle of virtual work applied to the body ji for a variation Su. in displacement, yields: 
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Figure 6.4 Separation of Cohesive zone from the bulk material. 
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[(j^,5u,dV = Q. 


( 6 . 8 ) 


divergence theorem to the above equation results in: 

\(T,jit^M')Ss,^{t + M)dV= \TXt + td)Su,{f + dJ)dS+ \T,{t+M'^,{t + td)dS, (6.9) 

where, r + Ar is the time of interest. Now defining the following terms: 


M^TXt+ht)-T,{tl 

(6.10) 

Au, =uXt+At)-uXt). 

(6.11) 


(6.12) 


(6.13) 


Substituting Equations 6.10-6.13 into Equation 6.9 results in the following: 


- I ^T,5^u,dS+\^a,J6^s,JdV= \T,{t-^i^)8iai,dS^ \ T,{f)5{Lu)dS- Ls,jdV . 
spi(i) p 3A a^ 2 (') p 


(6.14) 


Equation 6.14 represents the virtual work for a general body with initial and time dependant 
boundaries. Rather than applying this formulation directly to a finite element, a separation is 
considered (as shown in Figure 6.4) between the body y^and the boundary dP 2 ^\ snd the 
cohesive zone, with boundary dp^{t)- This separation yields: 



(6.15) 

1 

II 

(6.16) 

T,{np^{,))=<i"P2i>))- 

(6.17) 
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Thus, writing the incremental ized principle of virtual work for a cohesive zone results in: 

I LT,S^ix,dS+ \^a,^5Ks,^dV = - f T{t)5 tu,dS - Hs^dV . (6.18) 

Pc sp^it) P^ 

Assuming that the cohesive zone has no volume in the initial undeformed state and is completely 
described by boundary tractions, Equation 6.18 reduces to 

J AT,SAu,dS = - lT,{t)SAu,dS. (6.19) 

SPei') dp,(,) 

The traction-displacement relation may be written in the following incremental form: 

^T, = k,J^[tl,]+^T,\ ( 6 . 20 ) 


where the As represent changes in the associated quantities during the time step Af .The last term 
in the above equation in zero and the terms can be obtained by differentiating Equations 6. 1-6.3 

with respect to the displacement jump. 

Substituting Equation 6.20 into Equation 6.19 results in the following variational form: 


[ k,^liii)^Au,dS+ J AT,‘^SAu,dS = - 

dPcil) %(') ^P(‘) 

where, 



( 6 . 21 ) 

( 6 . 22 ) 

(6.23) 

(6.24) 
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In Older to simplify the integration, the surface area of the cohesive zone may be written in terms 
of an upper and lower surface: 


5y3c(/) = 5/?;(/)+s/?i(/). 


(6.25) 


Rewriting the principle of virtual work for the cohesive zone thus results in 




6.26) 


6.3.2 Finite element discretization 

Equation 6.26 may be discretized in space by writing the displacement field on each surface as a 
function of the nodal displacements on that surface. The crack opening displacement, [m, ] must be 
written with respect to each surface. This is derived from the definition of the traction- 
displacement relationship (in Mode I, [w,] for the lower surface will be negative). The nodal 

displacements on upper and lower surface are indicated as w“ (w), and m“ (/), resi>ectively. Thus, 
for the upper surface, it is assumed that: 


u, = N°u^(u), Au, = A^“Au“(m), (6-27) 

[«,]= , A[«,]=A^“Ai/r W-A^“Aw“(/). (6.28) 

Similarly, for the lower surface, it is assumed that: 

u, = Am, = N^Au^il), (6.29) 

[m,]=A^“m,“(/)-A^“m,“(«) . a[m,]=7/“Am,“(/)-A^“Am,“(m), (6.30) 
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where, N are the shape functions. Assuming a linear variation in the displacement field (a = 4) 
on the upper and lower surface, the interpolation functions for the mapped coordinate system are 
as given by Reddy [36]. A typical element is shown in Figure 6.5. Figure 6.6 illustrates the normal 
crack opening displacement as defined for the lower surface (with a linear variation in the 
displacement field). Substituting the interpolation functions into the principle of virtual work 
yields: 


Lo) W- M (/)]dZV“ Au“ {u)dS + {u)dS 


Because both the upper and lower surfaces are identical in the undeformed coordinate system the 
shape functions are equal. Equation 6.31 thus simplifies to the following: 


(»)- A»'(0K* KW-K + (“)<« 

+ (0- Au;(«)]w /iuf(!)ds + 


Factoring out the virtual displacements, SAuf («) and SAuf (/) , the principle of virtual work yields 
the following equilibrium equations for the upper and lower surfaces: 
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Figure 6.5 Parent element for interpolation and integration. 



Figure 6.6 Crack opening displacement in Mode I. 


90 



t"(,) dS=- dS , 

Ip/ (,)*„ [a«; (/)- Ai/P (m)J dS + Ip, ~ Ip/ (,)^(0 A^“ dS . 


(6.33) 

(6.34) 


The above may be written in matrix notation: 


k(»)]m'W}4„;(')))={F,M). (6.35) 

K(;)1{{a.,»(/))-4;(„))}={f,(;)}, (6.36) 


where, 






(6.37) 

(6.38) 


Rearranging and combining Equations 6.35 and 6.36 yields: 


■ k(')], 
[-^,(“)1 



[|^u5(m)Jj 1{a;(w)}J 


(6.39) 


Eq. (49) represents the local element sub-matrices in the p, n, and s directions. Expanding Equation 
6.39 in the p, n, and s directions one can obtain a general form: 


[^lower ] [“ Slower] | {^lo^er }1 _| {f^lower } 1 

["^upper] [A^upper] _ [{A//ipper}J \ Rapper }J 


(6.40) 
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where, and and J and represent the assembled stiffness matrices 

and force vectors for the lower and upper surface, respectively. In addition, and 

are the local displacement vectors. 


6.3.3 Global assembly 


One can expand Equation 6.40 in the p, «. and s directions to form element submatricies. For 
simplification in assembly, one only has to form the stiffness matrix associated with the lower 
surface. For example, if a = 4 

(Linear Interpolation), [.^y(/)j would be a 4 x 4 matrix; and a 12 x 12 matrix, would be 

assembled in the following manner: 

"a:),’ 0 0 0 0 a:;,^ 0 o a:],'* o o' 

0 ATi' 0 0 0 0 a:*^ 0 0 a:^'^ o 

0 0 a:]' 0 0 a:]2 0 o a:]^ o o a:]'* 

a:J' 0 0 Kf Q 0 Kf 0 0 Kf 0 0 

0 A:^' 0 0 ATf 0 0 Kf 0 0 Kf 0 

0 0 Kf 0 0 Kf 0 0 Kf 0 0 Kf 

Kf 0 0 a:J' 0 0 a:J' 0 0 ATf 0 0 

0 Kf 0 0 Kf 0 0 Kf 0 0 ATf 0 

0 0 a:,^' 0 0 Kf 0 0 Kf 0 0 Kf 

Kf 0 0 0 0 Kf 0 0 Kf 0 0 

0 Kf 0 0 Kf 0 0 Kf 0 0 Kf 0 

0 0 Kf 0 0 Kf 0 0 Kf 0 OKf 


Knowing [a:,„,.„] and assembling the force vector, one can write local equihbnum m 

terms of the lower surface: 
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required. Thus the main program uses the subroutine to construct the stiffhess matrices and load 
vectors of the cohesive elements in the model, and then assembles these matrices and load vectors, 
together with those of the bulk elements to construct the global equations. Incremental analysis 
with the Newton Raphson iterative scheme is used. 

Eight noded isoparametric element has been used as the cohesive element Since the upper 
and the lower surfaces of the element are co-planar in the initial configuration, and since the 
integration has to be performed only on the areas, four Gauss points are used. The stiffhess and 
internal load of the element is assumed to be contributed only by the Gauss points. Hence, the 
criterion of failure of the element is that, when X becomes unity at a particular Gauss point, the 
contributions of that Gauss point to the stiffiiess and load vector of the element vanish, and the 
element is assumed to have failed, when X exceeds unity at all the four Gauss points. 


6.4.1 One Element Tension Test Problem 

To test the subroutine, a simple tension test problem (as shown in Figure 6.7), with one bulk 
element (8-noded isoparametric) and one cohesive element, has been solved under displacement 
controlled loading. Displacement controlled loading is used because convergence problems arise 
with load controlled loading when the number of elements in the problem is less, due to the 
negative stiffness encountered during the softening of the cohesive elements. 

The nuclear reactor pressure vessel steel SA 508 has been used as the material. 
Deformation plasticity model has been used for the bulk elements, with Young’s Modulus equal to 
210 GPa, yield strength 450 MPa, hardening exponent 10.6, and yield offset equal to 0.91 [37]. 

The parameters of the cohesive element have been chosen in the following manner; is 

chosen same as the yield strength of the bulk material (450 MPa). The value of 5„ is chosen such 
that the initial stiffness of the cohesive element is comparable to that of the bulk material. This 
gives a value of 0.01 mm for 5„. Since the problem is purely of mode I, the shear terms have been 

ignored. 

The boundary conditions used are as shown in Figme 6.7. Based on Are above, the problen. 
is solved, and the traction.displacen,en. tesponse of the cohesive element, as obtained from Are 
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Normal Traction ^Pa) 


solution is ploltod along with the material constitutive curve in Figure 6.8. As expected, the 
cohesive element obeys the constitutive model used. 



Figure 6.7 One-element Tension test problem. 



Total Normal Separation (mm) ► 

Figure 6.8 Response of the cohesive element under Tension test problem. 


95 





6.4.2 Crack Growth Analysis of a DCB Specimen 

Figure 6.9 shows the Double Cantilever Beam (DCB) Specimen used for the analysis. The half 
thickness of the beam (h), as shown in the Figure is taken as 10 mm. The total length of the Beam 
is taken as 200 mm, and the width of the beam is taken equal to 100 mm. The value of the initial 
crack length (n,) has been taken to be 5 mm. The material properties and the values of the 
parameters of the cohesive zone are taken same as described in section 6.4.1. 

Only the symmetrical half of the specimen has been modeled (exploiting the symmetry 
about the crack plane), with a layer of 8-noded cohesive elements throughout the crack plane 
(except the initial de-bonded surface), as shown in Figure 6.9(b). The size of the cohesive elements 
in the direction of crack propagation is 1.67 mm. Thus, crack extension in increments of 1.67 mm 
is obtained. Also, 8-noded isoparametric elements have been used to discretize the bulk material. 
Again, a displacement controlled loading has been considered. The results of the analysis are 
presented in section 6.4.3. 



Figure 6.9 Crack growth analysis of a DCB Specimen. 
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6.4.3 Results and Discussion 


Figure 6.10 shows the crack propagation with increasing load-line displacement. The crack 
starts growing when the load-line displacement reaches a value of about 0.5 mm and continues to 
grow up to a crack length of about 29 mm, after which, the rate of crack growth with respect to the 
load line displacement becomes very small. Figure 6.11 shows the load-displacement relationship 
at the load line. Before the crack starts growing, the load increases linearly with the displacement. 
As the crack starts propagating, the load decreases due to the increase in the compliance with 
increasing crack length. 

Energy release rate (Gy) has been calculated using the load-displacement curve and is 
plotted in Figure 6.12. The nature of the curve matches very well with that of the crack-growth 
resistance curves available in the literature [22, 37]. 
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Figure 6.12 Crack growth resistance curve. 
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6.4.3 Crack Growth Analysis of a nozzle-corner crack 


The method is finally used for the crack growth analysis of a given no 2 zle-comer crack. A quarter- 
circular crack of 60 mm depth is taken as the initial crack configuration for analysis. The nozzle- 
vessel-intersection is meshed in ABAQUS in the same manner as described in Chapter 4, except 
that 8-noded isoparametric brick elements are used in the region near the crack, instead of 20- 
noded brick elements. In addition, 200 cohesive elements are incorporated in to the FE model near 
the crack, as shown in Figure 6.13. The region shown protruded in the Figure is additionally 
modeled here and first meshed with 8-noded isoparametric brick elements in compatible mode 
with the brick elements in the region near the crack. These elements are then declared as user 
elements in the input file generated by ABAQUS and the upper and lower surfeces of these 
elements are merged togedrer by editing the coordinates of the nodes on the upper (outer) surface. 

As seen in Figure 6.13, 40 cohesive elements are taken along the crack profile and 5 
elements are taken in the direction normal to the profile. The size of each cohesive element is 
about 2 mm along the profile and 1 mm normal to the profile. 


Vassal Side 



Nozzle Side 


Figure 6.13 Cohesive elements near the crack. 
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The material properties of the bulk material and the parameters of the cohesive zone are taken in 
the similar manner as described in section 6.4.1. Internal pressure load of 15 MPa is applied and 
the boundary conditions are taken same as described in Chapter 4. In addition, the upper surface of 
the cohesive zone is fixed in all directions. 

Load controlled incremental analysis is performed. The solution converged up to a pressure 
of 9. 1 MPa, and after that, the rate of convergence became very low, requiring very small time 
increments. The crack growth results obtained from the analysis are presented in the following 
section. 


6.4.4 Results and Discussion 

Figures 6.14-6.19 show the status of the cohesive elements (failed or surviving), at different values 
of internal pressure. The elements shown in red have failed, and the elements shown in green are 
still surviving. The elements in which all the four gauss points have failed, are considered to have 
failed. The boundary between the two colours thus represents the current crack profile. 

As seen in Figure 6.14, the crack has not yet started propagating at a pressure of 7 MPa. At 
a pressure of about 7.5 MPa, the crack has started propagating at the nozzle side. This is consistent 
with the SIF results obtained in Chapter 4, wherein it is seen that the SIF is higher at the nozzle- 
side of the crack. 

It is further observed form Figures 6.15-6.19 that the rate of crack propagation increases 
with the increasing pressure, once the crack has started growing. Finally, at a pressure of about 9.1 
MPa, the crack has grown by about 4 mm in the nozzle direction and about 1 mm in the vessel 
direction, thereby assuming almost a quarter-elliptical shape. This is again consistent with the SIF 
results, which showed that for a quarter-elliptical crack, the percentage variation of SIF is less than 
that for a quarter-circular crack, and hence an initially circular crack will grow in such a way to 
assume an elliptical shape. 

Since the total crack growth is of the same order of the size of the cohesive elements, a 
finer mesh can be expected to produce even better results. 
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Figure 6.14 Status of cohesive elements at internal pressure of 7 MPa. 
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Figure 6.15 Status of cohesive elements at internal pressure of 7.5 MPa. 
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Figure 6.1S Status of cohesive elements at internal pressure of 9 MPa. 
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Figure 6.19 Status of cohesive elements at internal pressure of 9.1 MPa. 
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6.5 Closure 


A three-dimensional formulation of the cohesive zone method for crack growth analysis has been 
presented. The method has been applied for the crack growth analysis of a DCB specimen and the 
results are found to match qualitatively well with the available literature. The methodology is then 
used for the estimation of stable crack growth in a given quarter-circular nozzle-comer crack under 
a constant internal pressure load (normal operational load). The pattern of crack growth has been 
found to be consistent with the pattern of variation of SIF along the profile of a quarter-circular 
crack as estimated in Chapter 4. The method can be conveniently extended to estimate the stable 
crack growth during a Pressurized thermal shock also. 
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Chapter 7 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 


7.1 Conclusions 

Structural Integrity of nuclear reactor pressure vessels with no 22 :le-comer cracks has been 
investigated under normal operational condition, and also under a loss of coolant accident. Linear 
Elastic Fracture Mechanics has been used for the purpose and the crack characterizing parameter. 
Stress Intensity Factor, has been evaluated along the profile of the nozzle-comer cracks for 
different geometrical and crack parameters. The results have been found to be within 4 % of 
accuracy. The individual effects of different parameters on Stress Intensity Factors have been 
studied. As an attempt to predict the stable crack growth of nozzle-comer cracks, the cohesive 
zone technique for three dimensional crack problems has been presented, and is applied for the 
crack growth analysis of a Double Cantilever Beam Specimen. The following are the main 
conclusions drawn from the present work: 

1 ) Use of quarter point elements developed by a simple shift of mid side nodes, rather than the 
face-collapsed wedge shape quarter point elements, to simulate the singularity near the 
crack, has been found to produce sufficiently accurate results. 

2) Under the normal operational load condition, sufficient factor of safety (greater than 5) 
exists, even for the half thickness cracks. Nozzle side of the comer-crack has a higher value 
of Stress Intensity Factor, and hence, there seems to be a possibility of Leak-before-Break. 
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3) The normalized average Stress Intensity Factor has been found to be almost independent of 
the crack aspect ratio. 

4) The Stress Intensity Factor has been found to bear simple relationships with all the 
parameters individually, thereby indicating the possibility of development of more accurate 
empirical relations for its calculation, than already available in the literature. 

5) Under the faulted condition, LOCA, the cracks sufficiently safe earlier, may become 
critical. Also, since the vessel side of the comer-crack has a higher value of Stress Intensity 
Factor during the faulted condition, the possibility of Leak-before-Break does not seem to 
be there. 

6) The cohesive zone technique has been realized to be a powerful tool for the estimation of 
stable crack growth. The consistency between the pattern of growth of the given nozzle- 
comer crack and that of SIF variation over the crack profile, suggests that the crack growth 
results are satisfactorily well. 


7.2 Scope for Future Work 

The present work provides some insight to the assessment of stmctural integrity of nuclear reactor 
pressure vessels having nozzle comer cracks, but the following work is also of interest and can be 
pursued as a future work: 

1) Elasto-Plastic Fracture Mechanics can be applied to capture the effect of plastic zone near 
the crack. J-integral can be evaluated over the crack profile, following an Elasto-Plastic 
stress analysis. 

2) Stress Intensity Factor and J-integral database can be generated, followed by 
development of empirical relations for their calculation. 

3) The possibility of crack arrest, in the event of a crack propagating after becoming critical, 

needs to be assessed. 

4) The amount of stable crack growth during a Pressurized thermal shock can be estimated, 
so that, the number of such transients that the pressure vessel can withstand before 
becoming operationally unsafe, can be predicted. 
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